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ABSTRACT 

We consider torsional oscillations of magnetars. This problem features rich dynamics 
due to the strong interaction between the normal modes of a magnetar's crust and 
a continuum of Magneto-Hydro-Dynamic (MHD) modes in its fluid core. We study 
the dynamics using a simple model of a magnetar possessing a uniform magnetic 
field and a thin spherical crust. Firstly, we show that global torsional modes only exist 
when one introduces unphysically large dissipative terms into the equations of motion; 
thus global modes are not helpful for understanding the magnetar Quasi-Periodic 
Oscillations (QPOs). Secondly, we solve the initial- value problem by simulating the 
sudden release of an initially strained crust and monitoring the subsequent crustal 
motion. We find that the crustal torsional modes quickly exchange their energy with 
the MHD continuum in the core, and decay by several orders of magnitude over the 
course of ~ 10 oscillation periods. After the initial rapid decay, the crustal motion is 
stabilized and several time- varying QPOs are observed. The dynamical spectrum of the 
simulated crustal motion is in qualitative agreement with that of the x-ray light-curve 
in the tail of a giant magnetar flare. The asymptotic frequencies of some of the QPOs 
are associated with the special spectral points — the turning points or edges — of the 
MHD continuum, and are not related to those of the crust. The observed steady low- 
frequency QPO at 18 Hz is almost certainly associated with the lowest frequency of the 
MHD continuum, or its first overtone. We also find that drifting QPOs get amplified 
when they come near the frequencies of the crustal modes. This explains why some 
of the observed QPOs have frequencies close to the expected crustal frequencies, and 
why these QPOs are highly variable with time. 



1 INTRODUCTION 



Magnetar oscillations have recently been the subject of 
intense theoretical interest. This interest is motivated by 
quasi-periodic x-ray luminosity oscillations (QPOs), which 
were observed in the tails of 2 giant Soft Gamma-ray Re- 
peaters' (SGR) flares (Israel et al. 2005, Strohmayer & Watts 
2005, 2006, Watts & Strohmayer 2006, see also Barat et 
al. 1983). The QPOs typically last for ~ 100 seconds, are 
detected with high signal-to-noise ratios, and have frequen- 
cies which range from 18 to 1800Hz. The QPOs open an 
exciting possibility to directly explore the physics of magne- 
tars, and their correct interpretation is of great importance. 

SGR QPOs have been commonly interpreted as pure 
crustal shear modes (Duncan 1998, Piro 2005, Lee 2006, 
Samuelsson & Andersson 2006, Watts & Reddy 2006, Sotani 
et al. 2006a). If this interpretation were correct, it would al- 
low one to measure or strongly constrain the shear modulus 
and depth of the crust, an unprecedented feat in the neutron- 
star astrophysics (Strohmayer & Watts 2006). However, the 
presence of the strong magnetic field which exists inside a 
magnetar may present difficulties for this interpretation. In 
particular, Levin (2006, L06) has pointed out that hydro- 
magnetic (MHD) mechanical coupling between the crust and 



the core occurs on the timescale < 0.1 seconds, and should 
be taken into account. L06 made 2 basic points: 1. MHD 
coupling ensures that pure crustal modes do not exist, and 
global modes of the whole star must be considered, and 2. 
Long-term survival of the global mode is in danger, since it 
is expected to couple to a continuum of MHD modes (the 
Alfven continuum) in the core, and this coupling would act 
to damp the mode. More recently, Glampedakis, Samuels- 
son, & Andersson 2006 (GSA) and Sotani et al. 2006b (S06b) 
have found global MHD-elastic modes in simple toy magne- 
tar models, and have argued that the analogues of these 
modes produce QPOs in real magnetars. However, both of 
the toy models have been constructed in a way which ex- 
plicitly excludes the presence of an Alfven continuum in the 



* GSA make use of the rectangular geometry with a uniform mag- 
netic field, which ensures that all of the Alfven modes with the 
same quantum number have the same frequency. S06b consider 
dipole magnetic field in the spherical magnetar, but explicitly ex- 
clude / ± 2 coupling in their equations. This enforces the spherical 
symmetry in the physical problem described by their equations 
of motion. In fact, their equations are identical to those that de- 
scribe oscillation of the star with a purely radial spherically sym- 
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The coupling of hydrodynamical waves to Alfven con- 
tinuum has been extensively studied in the context of so- 
lar corona, and is well understood (Ionson 1978, Hollweg 
1987a,b). The absorption of the waves by Alfven continuum 
is sometimes referred to as the resonant absorption. In this 
paper we build on the work done by the solar physics com- 
munity and undertake a thorough investigation of the influ- 
ence of the Alfven continuum on the oscillatory behavior of 
the magnetar crust. The plan of the paper is as follows: 

In the next section we describe our simplified, but topo- 
logically correct magnetar model. We derive equations of 
torsional motion, and search for normal modes of the sys- 
tem. We find the normal modes only exist if one introduces 
sufficiently large frictional forces, e.g. the ones of the form 
—'ydr/dt, into the equations of motion. The eigenfrequencies 
are complex, and the modes decay with the rate indepen- 
dent of 7, typically over several oscillation periods (i.e., a 
fraction of a second). However, the normal-mode analysis 
is inconclusive, since the real frictional forces in a magnetar 
may be small and the normal modes likely do not exist. Thus 
in section 3 we turn to the initial value problem. We model 
the continuum by 10 4 oscillators chosen to mimic closely 
the MHD dynamics of the core (one can think of this idea 
as similar to the one behind spectral codes in fluid dynam- 
ics, except that our equations are linear. The true linear 
behavior of the real magnetar is recovered when the number 
of oscillators tends to infinity). We begin the simulations 
by releasing an initially strained crust, and then monitor its 
motion. The result of one such simulation is shown in Fig. 9. 
We find that the crustal deformation energy is quickly con- 
verted into the energy of MHD continuum in the core, as 
was predicted in L06 and as is suggested by the large imag- 
inary components of the normal-mode frequencies obtained 
in section 2. The amplitude of crustal motion is reduced by 
10 2 over several oscillation periods, but is then stabilized as 
the crust reacts to the vigorously-moving core. In this sec- 
ond time-interval we find QPOs in the crustal mo- 
tion, see Fig. 10. The asymptotic frequencies of some of the 
QPOs are associated with the special spectral points of the 
continuum. Both turning points and edges of the continuum 
produce QPOs (see section 3 and Figure 4 for an explanation 
of what these are). We also find that when the frequency 
of a drifting QPO approaches that of a crustal mode, the 
QPOs amplitude gets significantly amplified. Thus crustal 
frequencies feature intermittently enhanced power, in agree- 
ment with the observations. 

In section 4 we present an outlook on the outstanding 
theoretical questions related to magnetar QPOs. 



2 BASIC MODEL AND ITS NORMAL MODES. 

Finding oscillatory modes of magnetic stars presents a 
formidable computational and conceptual problem. Rincon 
& Rieutord (2003) and Reese, Rincon, & Rieutord (2004) 
in a tour-de-force calculation have computed normal modes 
of an incompressible fluid shell threaded with a dipole mag- 
netic field. Partly motivated by their work, we choose a sim- 
ple magnetar model, with several basic assumptions: 1. We 
take the elastic crust to be thin compared to the core size. 
This is an excellent approximation for mode frequencies less 
than ~ 600Hz. 2. We assume that the fluid core has uniform 
density and is threaded by a homogeneous internal mag- 
netic field directed along the z-axis. This assumption does 
not change the physics of the problem, but does simplify the 
calculations and makes them more transparent. 

We also consider oscillations with purely asimuthal, 4>- 
independent displacements £ and £ of the core and the crust, 
respectively: 



£r = = £r = £ e 
^ = fr(r,0), 



o, 



(i) 



Here r, 6, <fr are the spherical coordinates. We have made use 
of the thin-crust assumption in writing the last equation. 

We are now in the position to write down equations of 
motion for small displacements of the crust and the core: 



d 2 



dt 2 



dt 2 



Lei (£<#>) + Lb, 

'd% 

dz 2 



dt 



(2) 
(3) 



where the partial derivative on the right-hand side of Eq. (3) 
is evaluated along the z-direction. Here c a is the Alfven ve- 
locity, L e i (£<#,) and Lb is the acceleration of the crust due to 
the elastic and magnetic stresses, respectively. We have in- 
troduced the frictional term —'yd^/dt into Eq. (3); we will 
show shortly that this term is needed to regularize the reso- 
nant response of the core to the periodic motion of the crust 
and is crucial for the existense of a normal mode. The ex- 
pression for Lb is derived below, while that for L e i is derived 
in the Appendix: 



Lel(^) = We 



oe 2 



+ cot(0)%-(cot 2 (0)-l)^ 



where u) e \ is the frequency given by 



R 



(4) 



(5) 



where p and p are the vertically averaged shear modulus 
and density of the crust, respectively, and R is the radius of 
the star. 



metric magnetic field, and thus, just as in GSA, the local Alfven 
waves with the same quantum numbers have the same frequency. 
We note that Rincon & Rieutord (2003) and Reese, Rincon, & 
Rieutord (2004) have previously found a continuum of torsional 
modes in the MHD configuration identical to that in S06b (they 
have not used any simplifying approximations in their analysis). 
To sum up, in both GSA and S06b the symmetry of their toy 
models collapses the Alfven continuum onto the discrete set of 
frequencies. 



2.1 Dynamics of the core: continuum of modes 
and response to periodic forcing. 

It is instructive to consider the motion of the core fluid, with 
the assumption of a fixed or periodically moving crust as 
an external boundary condition. The former will elucidate 
the structure of the Alfven continuum, while the latter is 
instrumental in the normal-mode analysis. 
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Hydro-magnetic stress enforces a no-slip boundary con- 
dition at the crust-core interface: 



(6) 



Let us consider the dynamics of the core, with the assump- 
tion of a fixed crust and zero friction. The core displays 
a continuum of singular oscillatory solutions to the pair of 
equations (3) and (6). These solutions are localized on cylin- 
ders of radius 770, with < r)o < R, and their mathemati- 
cal form is expressed most easily in cylindrical coordinates 
r/,z,4>: 



£<j,(r),z,<f>,t) = 5(r) - r]o)sm[nnz/h(r)o)] x 
exp[iimc a t/h(rio)] 



(7) 



£,4,(j),z,<j>,t) = 5(?7-r?o)cos[(n + l/2)7r«//i(?7o)] x 

exp[innc a t/h(rio)]; (8) 



cf. section 3 of L06. Here h(r)o) = yj R 2 — r) 2 , is the height 
of each cylinder, and n is an integer. While this continuum 
of MHD modes was derived for a simple magnetic-field ge- 
ometry, it must exist in other field geometries which can 
be obtained by continuous deformation of the uniform field. 
Moreover, since a continuos deformation of the field changes 
the mode frequencies continuously, the topology of the spec- 
trum remains unchanged. This means that the modes will 
form a countable set of one-dimensional continua; in other 
words, even for a complicated magnetic-field configuration 
the mode is parametrized by a pair of real and integer num- 
bers. The latter consideration will become important when 
we discuss general properties of QPOs in section 3. 

Now consider the core's response to a periodic motion 
of the crust, = g(6) exp(iwt), where lu could, in general, 
be complex. First we note that the geometry of our problem 
possesses the reflection symmetry with respect to the z = 
plane. Therefore the normal modes will be either even or 
odd with respect to z, and we restrict the crustal motion to 
that with g(9) = ±g(ir — 9). From Eqs. (3) and (6), we see 
that in the "odd" case the core motion is given in cyllindrical 
coordinates by 



while in the "even" case 



(9) 



(10) 



where 6(rj) = arcsin(?7/i?), h(rj) = y/ R 2 — rf , and k — 

y/iJ 2 - tUJ-y/Ca. 

2.2 Normal modes 

Now we are ready to derive the acceleration of the crust due 
to the hydromagnetic stress at the crust-core interface: 



Lb 



PcC a 



cos 9 



\dz), 



(11) 



where p c is the density of the fraction of the core fluid 
which participates in the Alfven motion, £ is the column 
density of the crust, and the partial derivative is evalu- 
ated at (rj,z) — (RsmO,Rcos9). By substituting Eq. (9) 



or (10) into Eq. (11), we get the following expressions for 
magnetically-driven acceleration of the crust: 

Lb{0) — — VaWi^r cot(wi cos6/va) cos0£<j,(6) (12) 



for the odd modes, and 

p R - 
Lb(0) = v a u)i-^- tan(u;i cos9/v a ) cos 9^(0) 



(13) 



for the even modes. Here v a — c a /R, and loi = uj 2 — iujy. 

We can now see that the Eq. (2), together with Eqs. (4) 
and (12/13), form an ordinary second-order differential 
equation for £$(9). The values of u> get selected by re- 
quiring that the solution of Eq. (2) satisfies the boundary 

conditions^ at the poles 9 — 0, 7r: 



zm = 



d9 2 



= 0. 



(14) 



We found it practical to solve Eq. (2) in the upper hemi- 
sphere, but require that at the equator either £$(9) = for 
the odd modes, or d^(9)/d9 — for the even modes. We 
also found it useful to make the substitution 



q{6) = im/o, 



(15) 



and rewrite the equations in terms of q(9). The new equa- 
tions do not have a singularity at the pole 9 — 0, and are 
very easy to integrate on the computer. We have checked the 
code by solving both analytically and numerically the case 
with Lb = 0. We find analytically that the wavefunctions of 
the free-crust vibrational modes are given by 



= dY l0 (9)/d9, 
with the eigenfrequencies 



(16) 



(17) 



LJi = - 1)(Z + 2)Wel. 

This /-scaling is in full agreement with that of Samuelsson & 
Andersson (2006). Our numerics gives excellent agreement 
with these results. 

Before we discuss our numerical results for the case of a 
non-zero Lb, one qualitative remark is due. From Eqs. (12) 
and (13), we see that if w is real and the friction coeffi- 
cient 7 = 0, then Lb diverges for the values of 9 which 
correspond to the location of the Alfven continuum mode 
in resonance with u>. It is these resonant interactions that 
are largely responsible for the exchange of energy between 
the global vibration and the Alfven cantinuum, and that 
are thus responsible for determining the imaginary part of 
uj (hence the name "resonant absorbtion" in the solar litera- 
ture). In our Runge-Kutta routine, we enforce small S-steps 
which scale as Lg 2 near the singularities (and we check that 
our results do not change when the step-size is reduced by 
a factor of 10). 

We now discuss our numerical results for normal modes 
of a magnetar. We find that the normal modes exist only 
when 7 is sufficiently large, i.e. 7 > 7 cr it- When this condi- 
tion is satisfied, then the frequency of the normal mode is 
complex and the mode decays with the rate close to 7 cr i t . By 



t These boundary conditions are derived mathematically from 
the Frobcnius series expansion of near the poles, or phys- 

ically by requiring the crustal angular velocity and acceleration 
be finite at the poles. 
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Figure 1. Damping rate for a global torsional mode is plotted 
as a function of the dissipative rate 7. The lines in the figure are 
numbered by an integer i = 1, 10; each lline corresponds to the 
magnetic field strength B c f{ = i X 10 14 G. For small 7, i.e. to the 
left of the lines, the modes don't exist. 



contrast, when 7 < 7 cr it, a thorough numerical search fails 
to identify the complex eigenfrequency for which all of the 
boundary conditions are satisfied. This is in full agreement 
with previous work done on the resonant absorption in the 
solar corona. For example, in Steinolfson's (1985) simula- 
tion the system behaves like a decaying normal mode when 
the friction is sufficiently large, while for small friction no 
normal-mode-like behaviour occurs and instead, the phase 
mixing is observed, where individual modes of the contin- 
uum are excited and oscillate each at their own frequency. 
The same occurs in our initial-value simulations, which we 
describe in Section 3. Hollweg (1987a,b) gives an excellent 
physical discussion of why, when the mode exists, its decay- 
rate is friction-independent^ 

In figure 1 we show the decay rate for the lowest- 
frequency normal mode as a function of 7, for the ten values 
magnetic-field strengths B c g = ^/A-k pc a ranging from 10 14 
tp 10 15 Gauss. 



t Briefly, this can be understood as follows: the energy is ab- 
sobed in narrow resonant layers; in our case they have cylindri- 
cal shapes. Friction produces 2 main effects: (a) it reduces the 
excitation amplitude of the resonant layer, and (b) it increases 
the effective width of the layer. It is easy to check that for the 
simple frictional terms we are using, the 2 effects compensate 
each other and the total absoption power is friction-independent. 
Hollweg (1987) proves that the same result holds for more com- 
plicated dissipative effects, like viscosity or Ohmic dissipation in 
the plasma. 



Here B c ff = \JBB C if the protons form a superconduc- 
tor with the critical field strength B c , and B c ff if the protons 
form a normal fermi liquid, and p is the density of the core 
material coupled to the magnetic field§. Our fiducial pa- 
rameters were pR/T, — 10 [probably an underestimate, but 
its increase would only increase the crust-core coupling — see 
Eqs (12) and (13)], uj e \ = 2n x 20Hz, and z/ a /o> e i = 0.2, ap- 
propriate for p ~ 10 14 g/cm 3 . We see that when the modes 
exist they show rapid decay, on the timescale <C 1 second. 
This result holds for all higher-order modes we have consid- 
ered, and shows how efficiently the crustal motion is coupled 
to the continuum. However, the friction (e.g., due to viscos- 
ity) may be small in magnetars, and the normal modes most 
likely do not exist. Thus the problem of magnetar torsional 
motion must be addressed using initial-value simulations. 
This is the subject of the next section. 



3 INITIAL- VALUE SIMULATIONS 

The aim of this section is to simulate torsional motion of a 
magnetar. During this motion, the discrete torsional modes 
of the crust interact strongly with a continuum of Alfven 
modes in the core, and this interaction affects dramatically 
the motion of the crust. In the next subsection, we explore 
with a help of a toy model the dynamics of a harmonic oscil- 
lator coupled to to a continuum of oscillators. The toy model 
provides us with an insight into QPOs of such a system, and 
gives us intuition for what to expect in the case of a magne- 
tar. In subsection 3.2 we set up the initial-value simulation 
for our magnetar model (the "real" magnetar, as opposed 
to the toy model in subsection 3.1), and present results. 

3.1 Coupling of a harmonic oscillator to a 
continuum: toy model. 

In Figure 2 we show the set-up of our toy problem. We 
consider a pendulum weighing 1kg, with a proper oscillation 
angular frequency of u)o = 1 (the units are irrelevant) . We 
model the pendulum coupling to a continuum of modes by 
suspending 10000 tiny pendulae, each weighing O.Olg, from 
the big pendulum. We arrange the angular frequency of the 
small pendulae to be 

to m = 0.5 + .0001 * m, (18) 

where m = 1, 2, 10000. Thus the frequency of the big pen- 
dulum lies in the middle of the range of those of the small 

§ This density can range from the density of protons p pr if the 
neutron superfluid is entirely decoupled from the proton motion, 
to the full core density p CO rc if the neutrons are efficiently cou- 
pled to protons. Our guess is that in reality p takes the value 
somewhere in between the 2 extremes: it may be hard for neu- 
tron superfluid to become completely decoupled from the charged 
component, since neutron superfluid vortices are expected to bc 
strongly magnetized and may interact strongly with the supercon- 
ducting flux-tubes; see, e.g., Ruderman et al. 1998 and references 
therein. We take the numerical value p = 10 14 g/cm 3 , about 1/10 
of the core density and twice the proton density. Our choices for 
B and p affect directly the value of the lowest-frequency QPO 
simulated in section 3. 
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Figure 2. Big pendulum coupled to a large number of small one. 
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Figure 3. Big-pendulum displacement as a function of time. 



pendulae. The initial condition for our simulation is as fol- 
lows: the big pendulum is deflected by a small angle (we 
keep the problem linear), while the small pendulae are re- 
laxed and hanging straight down. This is meant to mock an 
initially strained crust and relaxed core. The big pendulum 
is released, and the evolution is followed by two independent 
numerical routines. One routine uses the 4-th order Runge- 
Kutta method, while the other one uses a symplectic second- 
order leapfrog algorithm, which is very robust for simulating 
Hamiltonian systems (see, e.g., Kinoshita, Yoshida, & Na- 
gai 1991). Both runs conserve the total energy of the system 
with high precision, and produce results which are in ex- 
cellent agreement with each other. In Fig. 3 we plot the 



Figure 4. Zoom-in on the post-decay part of Fig. 3. Quasi- 
periodicity is clearly visible. 



big-pendulum displacement as a function of time. After sev- 
eral oscillations, the amplitude of the pendulum's motion is 
reduced by ~ 100, as the energy is rapidly transfered from 
the big pendulum to the small ones. Then the exponential 
decay abruptly stops, as the big pendulum now reacts to 
the collective pull of the small ones. The blow-up of this 
second region is shown in Fig. 4. The amplitude still decays, 
but only slowly, as 1/t. Astonishingly, even with naked eye, 
one can detect QPO(s) in the pendulum motion! Figure 5 
shows the time Fourier transform of the big-pendulum dis- 
placement for this interval of time. Two QPO frequencies are 
clearly present, 1.5 and 0.5, both identified with the edges 
of the continuum and NOT with the natural frequency of 
the big pendulum! 

We get a clue for the origin of these QPOs by plotting 
the phases of small pendulae as a function of the oscillator 
frequency, u) m ; see Figure 6. Over the range of m, the phases 
average out, thus preventing the small pendulae from pulling 
coherently on the big pendulum. The only location where 
this averaging does not occur is near the end points of the 
spectrum; see the arrows on Figure 6. Thus the pendulae 
near the end points of the spectrum do pull coherently on 
the big pendulum and produce the 2 QPOs observed in the 
simulations. The number of "coherent" oscillators shrinks as 
1/t, as the phase gradients with respect to m grow linearly 
with time. This explains why the QPO amplitude decays as 
1/t. In Appendix B we present a more mathematical way to 
understand this phenomenon. 

Apart from edges, there may be other special points in 
the continuum which could generate QPOs. One example 
is the local maximum or minimum in uj m as a function of 
m; we shall call such special points the turning points. The 
same reasoning as that for the edges, shows that the phases 
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Figure 5. Power spectrum of the post-decay big-pendulum dis- 
placement. Two QPOs are clearly visible; they are associated with 
the edges of the continuum frequency interval. 
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Figure 6. Phases of the small pendulae. 

of small oscillators near the turning point will not average 
out for some time, and hence these oscillators will act coher- 
ently. Moreover, the density of states is higher near a turning 
point than that near an edge, and diverges as \u — wo| -1 / 2 , 
where loo is the frequency of the reversal. We thus expect 
that turning points generate stronger QPOs than the edges; 
this is confirmed by our simulations shown below. The simu- 
lations also show that the turning-point-generated QPOs are 
longer lived than the edge-generated, and their amplitudes 
decay only as £ -1 / 2 . This decay law is explained mathemati- 



Figure 7. Spectral law with the turning point. 

cally in Appendix B. Figure 7 shows an example of a spectral 
law with a turning point. 

We simulate the initial-value problem for this example 
and show in Fig. 8 the big pendulum's displacement as a 
function of time. The strong long-lived QPO at the turning- 
point frequency of 0.5 is apparent with the naked eye. 

3.1.0.1 Effect of viscosity. We have modeled the effect 
of viscosity by introducing frictional forces between neigh- 
bouring small oscillators. We expect that with the passage 
of time the small oscillators get more out of phase, the ve- 
locity shear increases and so does the dissipation rate. This 
is confirmed by our simulations: the total mechanical en- 
ergy of the system is drained efficiently after some time. In 
magnetars, this may efficiently heat the core and affect the 
post-burst afterglow; we shall discuss these issues elsewhere. 
We find, though, that the QPOs survive for a much longer 
time than the global mechanical energy, because the oscil- 
lators creating it are precisely those ones which are moving 
in concert. The QPOs generated by the turning points are 
particularly robust. 

3.2 Initial-value problem for the magnetar model 

Lessons learned from the toy model lead us to expect (a) 
rapid decay of initial crustal perturbation and excitation of 
the core continuum, and (b) QPOs generated by the edges 
and turning points of this continuum. In our model for the 
magnetar we have turning points in all Alfven overtones, at 
angular frequencies tu n — mvc a /R for the odd modes, and 
u„ = (n — l/2)nc a /R for the even ones; the odd torsional 
motions are decoupled from the even ones. Thus potentially 
we expect QPOs at all of these frequencies. We shall find 
however, that our magnetar model displays a much richer 



© 0000 RAS, MNRAS 000, 000-000 



Magnetar QPOs 7 




Figure 8. Big-pendulum displacement as a function of time. Here 
it is calculated for the case when the small pendulae follow the 
spectral law in Fig. 7. QPO associated with the turning point is 
clearly visible in the post-decay time interval. 



dynamics than the toy models of the previous subsection, 
although the basic features of the toy models (initial rapid 
decay of the crustal motion, QPOs associated with contin- 
uum turning points) remain. In what follows we explain how 
our initial- value simulations are set up and show the results. 

3.2.0.2 Modal decomposition. A crustal displace- 
ment could be represented as a sum of the crustal normal 
modes: 



&(0,t) = X j b j (t)fM, 



(19) 



where fj(0) are proportional to the functions given in 
Eq. (16), and are normalized so that 



f 71 
Jo 



Me)f 3 (8)de = s K 



(20) 



We now develop a formalism which allows us to numeri- 
cally compute the time-evolution of bi(t). When the crust 
is not coupled to the core, b»(£) oscillate harmonically with 
the fruequencies of corresponding crustal modes, but their 
behaviour is very different when the crustal modes are cou- 
pled to the continuum modes in the core. In what follows 
we model this continuum with a large but finite number of 
small oscillators, and we test that our results do not change 
when the number of oscillators is increased. 

Consider, for concreteness, only odd modes (they re- 
main decoupled from the even ones because of the reflection 
symmetry). Recall that for a fixed cyllindrical radius r/o, the 
equation of motion for the core fluid is 



d ^(r)o,z) _ 2 d £</>{r)o,z) 

Co. 



dt 2 



dz 2 



(21) 



with boundary conditions £<t,(rio,h) = £?i[#(7?o)] and 
^4>(Vo, 0) = 0, where h — R 2 — rfi is the cyllinder's height, 
and 0(r)o) = arcsin(77o/-R). Let us introduce a new variable, 



x(Vo,z) = U(Vo,z) - ^4,[9{rj Q )]-. 

We obtain an inhomogeneous equation for x, 

d 2 x(yo,z) 2 d 2 x(Vo,z) _ z d 2 ^[e(- no ) 

dt 2 Ca dz 2 h dt 2 

but with easy boundary conditions 

x(Vo,0) = x(Vo,h) = 0. 



(22) 



(23) 



(24) 



Because of this boundary conditions, we can expand x m a 
Fourier series: 

x(vo, z, t) = £^° =1 a n (7/ , t) sin(n7rz), (25) 

where z = z/h(rj ). The right-hand side of Eq. (23) can be 
expanded in the same Fourier basis using the identity 

(-1)" 



z = -2£° 



■ sin(n7rz) 



(26) 



Now, by substituting Eqs. (25) and (26) into Eq. (23), we 
obtain the equations of motion for a n : 



d 2 a n (r]o,t) 
dt 2 



+ n 2 -K 2 v 2 a a n (r) ,t) 



2(~l) n d 2 U%o)] ^ 
im dt 2 



(27) 



We now use Eq. (11) to write down the expression for the 
hydromagnetic backreaction on the crust: 

L B {0) = -?§(ca/R) 2 (E~iror(-l)"a„fo(0)] +£(*)) ■ (28) 

The crustal mode amplitudes b m (t) obey the following equa- 
tions of motion: 



d 2 b„, 
dt 2 



+ ^mbm 



/>7T 

.A) 



L B (0)fm(0) sin 6d6, 



(29) 



where LO m is the frequency of the m's crustal mode. In writ- 
ing the last equation, we have used the normalization prop- 
erty of the modal wavefunctions f m {6). So far we have not 
used any approximations. Now we discetize the intergral 
in Eq. (29) by summing over a large number iV of points 
6i, which we take to be equally spaced with the interval 
A9 = tt/N. When we do the sum, the continuum mode apli- 
tude a n (8) is substituted with the discrete one, a; n = a„(0i). 
Thus effectively doing the sum instead of the integral substi- 
tutes a continuum of the core modes with the large number 
of the discrete core modes. The true continuum dynamics is 
fully recovered when N goes to infinity. 

We are now in the position to write down the equations 
of motion for the coupled crustal and core modes. From 
Eq. (27), we have 

d 2 a,i. 



+ n 7T v a a,i 



2(-l)\ 



dt 2 nn 
Further, from Eq. (29) we get 

d 2 b m ( 2 pR 

dt 2 



d 2 bm{t) 

= 1 Jm(Ut)- 



(30) 



f 2 . PR 2 
I OJm + —V a 



-^i^AflE„,iTMr(-l)" / m (0;) sin 6Ui. 



(31) 



Equations (30) and (31) are the main equations of this sec- 
tion. As with our toy models, we integrate these equations 
using 2 independent numerical techniques, the fourth-order 
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Figure 9. Crust displacement as a function of time. 



Runge-Kutta and the second-order leapfrog. The 2 meth- 
ods give results which are in excellent agreement with each 
other. 

We can now show some results from our numerical ex- 
periments. Here, we consider the 2 lowest odd crustal modes, 
bi and b 2 , with frequencies of 40Hz and 84.5Hz, coupled to 
10000 odd modes of the continuum: 1000 values of 9i with 
10 Alfven overtones at each point. This model should be 
representative for variability below 100Hz. For higher fre- 
quencies it is desirable to move away from the thin-crust 
approximation; this is the subject of future investigations. 
We begin with the crustal displacement bi = bi = 1 and 
the relaxed core, a;„ = 0. We monitor the crustal displace- 
ment at 6 = 7r/4. In Fig. 9 we plot this displacement as a 
function of time for the first few seconds after release. Like 
in our toy models, we observe an initially rapid decay of 
crustal motion, due to pumping of the crustal energy into 
the core. The crustal motion is then stabilized as the crust 
reacts to the vigorous movement of the neutron-star core. 
In Fig. 10 we plot the dynamical spectrum of the crustal 
motion for the first 100 seconds (we split the time-axis into 
1/2-second intervals, and take a Fourier transform at each 
interval. The density of points represents the magnitude of 
the square of the Fourier transform). In the dynamical spec- 
trum, we see several QPO-type features. The low-frequency 
QPOs asymptote to the turning points of the core contin- 
uum, c a nn/R (here we are considering the odd modes). In- 
termittent drifting QPOs appear at higher frequencies, and 
they get strongly amplified near the crustal frequencies. The 
nature of the QPO frequency drifts is unclear to us at this 
point. 

The simulated dynamical spectrum of Fig. 10 is in qual- 
itative agreement with that of the x-ray lightcurve in the tail 
of a giant flare, see Israel et al. 2005. In both cases there is 



Figure 10. Dynamical spectrum of the crustal displacement. 
Thin horizontal lines mark the pure crustal frequencies. Two low- 
frequecy QPOs asymptote to the continuum turning points, and 
are unrelated to the crustal modes. Larger version of the figure, 
with better resolution near the crustal frequencies, is available 
upon request. 



significant steady power below the lowest crustal frequency, 
and we conclude that the observed 18Hz QPO is almost cer- 
tainly the turning point of the core continuum. We also see 
an intermittent excess of power near the crustal frequencies, 
in qualitative agreement with the observations. 



4 OUTLOOK 

In this paper we have elucidated the crucial role that the 
core Alfven continuum plays in the dynamics of torsional 
motion of magnetars. We have shown that the global tor- 
sional modes do not exist unless the friction is unphysically 
large. We have performed a series of initial-value simulations 
for a simple but geometrically realistic magnetar model, and 
have observed QPOs with the properties closely resembling 
those in the tails of giant magnetar flares. In our model, the 
steady low-frequency QPOs are associated with the turning 
points of the Alfven continuum. This gives us a constraint 
on the combination of magnetic field geometry and strength 
and density of the core material which is coupled to the 
magnetic field. For our geometry, 



(io 15 g) 



10 14 g/cm 3 



1/2 



(32) 



Some of the higher-frequency power is clustered around the 
crustal frequencies, in agreement with the observations. This 
seems to be due to intermittent amplification of the drifting 
QPOs when they are close to a crustal frequency. 
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While the qualitative agreement between our simula- 
tions and observations is good, we see several directions for 
future research: 

1. Investigate qualitatively the origin of the QPO drifts is 
the simulations, and the reason why the QPO amplitudes 
get amplified near the crustal frequencies. 

2. Study the Alfven continuum for more realistic field geome- 
tries, e.g. the ones proposed in Braithwaite & Spruit 2004, 
2006. 

3. Investigate quantitatively the effect of viscous friction in 
the core. We have done some pilot studies for the toy models 
in the subsection 3.1, and have found that viscous friction is 
very efficient in draining of the core's kinetic energy, but does 
not significantly affect the QPOs. It is interesting to learn 
whether some of the long-term thermal afterglow could be 
generated from the heat deposited in the core by viscously 
damped Alfven waves. 

4. Investigate the dynamics of the magnetosphere. It is 
likely that the magnetosphere also features the continuum 
of Alfven modes, and they will affect the emission of x-rays 
in the giant flare afterglow. 
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6 APPENDIX A: ELASTIC FORCES AND 
ACCELERATIONS 

The purpose of this Appendix is to derive the expression (4) 
in the text for the crustal acceleration due to elastic forces. 
Since our problem is axially symmetric, the quickest deriva- 
tion is obtained by considering the flow of the z-component 
of the shell's angular momentum density. The displacement 
£</>(8) causes the horizontal shear stress in the shell which is 
given by: 



T 6i = WR)smO 



3(&(0)/sin0) 



(33) 



where /i is the shear modulus of the shell. The flow of angular 
momentum in the de direction is given by 



Lflux = 27r(i?sin(9) 2 / drT f 



(34) 



where the integration is over the thickness of the shell. The 
small thickness of the shell warrants the assumption of r- 
independent £$(0), so only /j, needs to be integrated. The 
angular momentum density with respect to the R9 coordi- 
nate is given by 



Ldensity = 2n(Rsm8) 2 T,(d£ lt> (6) /dt). 
Angular momentum conservation demands 
d(Ldensity)/<9i = -(l/.R)0(Lflux)/00. 



(35) 



(36) 



By substituting Eqs (34) and (35) into Eq. (36) and per- 
forming some trivial algebra, one obtains 



where 



de 2 



+ cot (61)- 



(cot 2 (60 - 1)& 



2 _ f v dr _ /Vp 



, (37) 



(38) 



The right-hand side of Eq. (37) is the operator L(£$) used 
in the text. 



7 APPENDIX B: QPOS FROM EDGES AND 
TURNING POINTS IN THE CONTINUUM 

The purpose of this Appendix is to explain mathematically 
the late-time behaviour of QPOs produced by the edges and 
turning points of the continuum in the toy models of sub- 
section 3.1. 



7.1 Edges 

Lets assume that after initial excitation, the oscillation am- 
plitude of the small pendulae described in the subsection 3.1 
does not change. Then the force acting on the big pendulum 
can be written as 



F(t) 



-L 



d,Lof(co) exp(iu)t), 



(39) 



where ^min and the upper and lower edges of the 

continuum, and function /(w) encompasses the amplitude of 
excitation, the coupling strength, and the density of the con- 
tinuum states at frequency ui. Lets continue f(ui) smoothly 
outside the (w m in, w m ax) region in such a way that / — > suf- 
ficiently fast (to be presently specified) as |w| — > co. We call 
this new function f(uj) and choose it in such a way that its 
inverse Fourier transform, f*(t), has a finite effective width 
At. Then F(t) is the convolution of /*(£) and s*(t), 



F(t)oc / drr(r)s*(t-T). 
Here 

s*(t) = [exp(iu; max t) - cxp(iaj m i n t)]/t 



(40) 



(41) 



is the inverse Fourier transform of the function s(w), which 
equals 1 for ii) m i„ < lo < a; max , and otherwise. When t 3> 
At, we can substitute l/(t — r) with 1/t in Eq. (40). Then 
we have for late times 

F(t) oc [/(Umax) cxp(icj max t) — /(u; m in) exp(iw m i n t)]/t , (42) 
and thus the edge QPOs decay as 1/t. 



7.2 Turning points 

Consider a downward turning point at angular frequency loo 
(i.e., a local maximum in a; as a function of m). The density 
of states near the turning point scales as (uio ~uj)~ l ^ 2 . Then 
the contribution of the oscillators near the turning point 
frequency to the force acting on the big pendulum is given 



by 

F{t) oc 



r 



duj(uiQ — lu) 1 ' 2 exp(iwt) 



(43) 
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Introducing the new variable x = (u)o — of)t and noticing 
that for large t the range of x becomes essentially (— oo,0), 
we see that 



Thus, the amplitude of the QPO associated with the turning 
point decays as t^ 1 ^ 2 . 
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